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1 .  Introduction 


Trajectory  and  advection  models  (Muench  1983;  Muench  and  Chisholm  1985; 
Muench  1989)  and  automated  display  systems  such  as  McIDAS  and  PROFS  (Schlatter 
et  al.  1985)  show  promise  of  providing  added  guidance  for  the  very  short  term 
forecasting  range  of  0  to  6  hours.  Full  synoptic  models  such  as  the  LFM  and  NGM 
have  skill  from  about  12  hours  to  beyond  48  hours.  However,  both  of  these 
approaches  fall  short  in  the  6  to  18  hour  period  which  is  perhaps  the  most  vital  to 
the  terminal  forecaster.  It  is  in  this  period  that  meso-/3  scale  disturbances  (such  as 
fronts  or  squall  lines)  which  are  in  the  local  area  will  have  an  immediate  impact  on 
the  terminal  forecast.  Advection  models  begin  failing  after  a  few  hours  because  of 
synoptic  and  mesoscale  changes  in  wind  patterns  and  because  topographic  and 
geographic  influences  are  often  not  incorporated  (Muench  1983;  Muench  1989).  On 
the  other  hand,  the  initialization  of  synoptic  scale  models  is  based  on  synoptic  scale 
data  which  does  not  normally  include  mesoscale  features  present  at  the  time  of 
initialization.  Hence,  the  synoptic  scale  model  cannot  be  expected  to  forecast  these 
mesoscale  features. 

Several  research  mesoscale  models  exist  which  are  capable  of  providing 
forecasts  on  the  desired  space  and  time  scales  (see  Pielke  1984,  Appendix  B,  for  a 
partial  list  of  current  research  mesoscale  models).  However,  these  models  tend  to 
be  just  as  computationally  intensive  as  operational  synoptic  scale  models,  which 
limits  their  practicability  in  an  operational  setting  —  each  region  of  the  country 
would  require  a  dedicated  supercomputer  to  provide  mesoscale  guidance  for  that 
region.  The  research  reported  here  seeks  to  provide  regional  mesoscale  guidance 
using  a  different  approach  which  is  well-suited  for  an  operational  setting.  A  meso-0 
model  is  being  developed  which  is  capable  of  providi'ig  an  operationally  useful 
forecast  on  a  super-micro  computer  (Vax  3  class).  Such  a  model  could  then  be  run 
"in-house”  at  each  forecast  office  to  provide  guidance  for  a  region  centered  on  that 
office.  This  approach  requires  the  formulation  of  a  model  which  is  considerably 
different  from  current  research  mesoscale  models.  The  goal  is  to  develop  a  model 
which  has  many  fewer  levels  and  less  complicated  physical  parametenzations  than 
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research  models,  but  which  is  stilt  capable  of  reproducing  the  physical  processes 
important  to  mesoscale  systems. 

The  next  section  will  describe  the  mathematical  formulation  of  the  model 
being  developed  here  including  the  numerical  treatment  and  physical 
parameterizations.  Section  3  will  present  results  of  tests  which  have  been 
conducted  with  various  prototype  versions  of  the  model.  The  radiation  and 
boundary  layer  parameterizations,  which  are  crucial  to  the  unique  aspects  of  the 
model  being  developed,  are  described  in  some  detail  in  section  4.  Section  5  gives  our 
conclusions  at  the  end  of  the  first  year  of  model  development  and  outlines  the 
research  which  will  be  continuing  in  the  future. 

2.  Model  Description 

a.  Basic  Model  Equations 

Mesoscale  models  capable  of  simulating  atmospheric  phenomena  have  existed 
for  more  than  a  decade  (Anthes  and  Warner  1978;  Nickerson  1979),  and  the  equation 
set  in  cr-coordinates  [cr==ip  pt  ps  — pt)  where  p,  is  the  surface  pressure  and  p,  is 
a  pressure  level  specified  as  the  top  of  the  model  (we  take  p,  =  100  mb  for  this 
study)!  is  well  established.  As  discussed  above,  however,  these  models  would 
require  a  supercomputer  environment  in  order  to  run  operationally.  To  develop  a 
mesoscale  model  which  can  run  operationally  on  a  relatively  small  computer,  we 
must  take  a  substantially  different  approach. 

The  two  obvious  ways  in  which  to  increase  the  speed  of  a  three-dimensional 
model  reducing  the  number  of  layers  or  decreasing  the  domain  size  —  are 
clearly  inadequate  in  and  of  themselves.  The  first,  by  itself,  would  lead  to 
inadequate  treatment  of  processes  near  the  surface  because  the  boundary  layer 
would  no  longer  be  resolved.  The  second  would  decrease  the  predictive  time  for 
any  mesoscale  phenomena  not  produced  by  local  forcing  since  it  is  impossible  for 
mesoscale  disturbances  to  advect  into  the  domain.  We  do  seek  to  decrease  the 
number  of  model  layers  in  order  to  decrease  the  number  of  computations  per  time 
step,  but  we  avoid  problems  near  the  surface  by  treating  the  boundary  layer 
explicitly. 
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Much  of  the  local  forcing  for  mesoscale  processes  takes  place  in  the 
planetary  boundary  layer.  Previous  mesoscale  models  packed  many  layers  near  the 
surface  to  resolve  boundary-layer  thickness  and  to  adequately  calculate  fluxes 
there  (see,  for  example,  the  ^-coordinate  model  of  Nickerson  1970).  We  have 
developed  a  two-layer  model  whose  lowest  layer  represents  the  boundary  layer. 
This  layer  is  not  fixed  in  pressure  or  a  and  is  allowed  to  change  its  depth  as 
function  of  time  representing  the  physical  structure  of  the  variable  depth 
boundary  layer.  The  fluxes  from  the  surface  into  this  layer  as  well  as  the  fluxes 
from  the  boundary  layer  into  the  layer  above  are  calculated  directly.  A  new 
concept  here  is  the  two-way  interactive  nesting  of  this  two-layer  model  into  the 
lowest  layer  of  a  four-layer  a- coordinate  model. 

The  equation  set  for  the  four-layer  mode!  is  essentially  identical  to  that  of 
Anthes  and  Warner  (1978).  1  he  map  factor  has  not  been  included  in  our  set  because 

a  scale  analysis  indicates  this  term  is  not  significant  for  the  meso-0  scale  domains 
we  will  be  considering.  The  lowest  layer  of  the  four-layer  model  is  on  the  order  of 
200  mb  thick.  It  is  defined  by  the  surface  (cr  —  1 )  and  the  constant  a  surface  at 
cr— cr*.,,  (see  F;ig.  1).  The  two-layer  model  is  nested  within  this  layer,  with  the  oh 
surface  which  divides  the  two-layer  model  representing  the  variable  depth  of  the 
planetary  boundary  layer  as  a  function  of  r,  y,  and  t.  Since  <7^  is  variable  in  sigma, 
computations  in  the  two-layer  model  require  a  coordinate  transformation  into  a 
system  in  which  a-  represents  a  constant  coordinate  surface.  We  call  this  new 
coordinate  system  “boundary  layer  coordinates"  and  denote  it  with  the  variable  7. 
We  define  7  as 


n  =  CT-_<7> 

H  ■ 

|crh 

-C utm 

a  ■ 

C . 

H 

l  1 

i  Vh 

cr 

(2.1) 


Prognostic  equations  for  a  and  i>  momentum,  t  (=-p*  p.),  and  temperature,  as  well 

as  diagnostic  equations  for  vertical  velocity  in  the  7  system  and  the  hydrostatic 
relation  become: 
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where 

V  =  vertical  velocity  in  ^-coordinates 
u>  —  pressure  vertical  velocity  dp  dt ) 

<7ibrri  =  vertical  velocity  at  the  dkt;  interface 
=  geopotential 
F  =  friction  term 

R  =  specific  gas  constant  for  dry  air 

cp  —  specific  heat  of  dry  air  at  constant  pressure 

/  =  Coriolis  parameter. 
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In  (2.4i  and  (2.6),  H  and  H refer  to  the  values  of  H  in  the  upper  and  lower  layers 
of  the  two-layer  model,  respectively  (see  (2.1)1. 

The  pressure  vertical  velocity,  u,  is  found  through  a  diagnostic  relation  for 
each  V  layer 


u j  =  < i? H  ‘-cTn’ 


dt 


4; 


•T 


-  [T1  rl) 


-ct2» 

3t 


r  u 


da b 
3r 


da^ 


t  rjH 


(2.8) 


where  the  f  or  sign  in  the  second  term  on  the  RHS  is  chosen  for  the  upper  and 
lower  layers,  respectively.  Note  that  this  expression  includes  explicitly  the 
temporal  and  spatial  variation  of  the  boundary  layer  height,  <7n.  The  boundary 
layer  height  is  a  product  of  the  boundary  layer  parameterization  (see  section  4), 
and  the  rate  of  change  of  this  height  Ocr^  ?t )  is  also  computed  in  the 
parameterization  scheme  for  use  in  (2.8). 

While  cr-coordinate  prognostic  variables  are  weighted  by  t  in  the  typical 
formulation  (see  Anthes  and  Warner  1978),  it  is  obvious  in  the  above  equations  that 
^-coordinate  variables  are  weighted  by  ^ H .  Other  prognostic  equations  similar  to 
(2.5)  can  be  written  to  provide  tendencies  for  trace  variables  such  as  specific 
humidity. 


6.  Vertical  and  horizontal  nesting  in  the  model 

The  merger  of  the  two-layer  model  into  the  four-layer  model  represents  a 
classical  two-way  interactive  nesting  problem.  The  two-layer  model  actually 
replaces  the  lowest  layer  of  the  four-layer  model.  Information  from  the  upper 
three  layers  of  the  four-layer  model  is  passed  to  the  two-layer  model  through 
vertical  differences  across  the  interface  separating  the  models  as  well  as  through 
the  downward  integration  of  divergence  in  the  calculation  of  change  in  surface 
pressure  (the  first  term  on  the  RUN  of  (2.4)].  The  four-layer  model  also  provides 
o,.bm  for  use  in  (2.6).  Information  is  passed  back  from  the  two-layer  to  the  four- 
layer  model  through  vertical  differences  across  the  interface  and  through  the  T- 
weighting  of  variables  in  the  four-layer  cr-coordinate  model  (since  T  includes 
contributions  from  all  layers).  Kvcn  though  no  prognostic  equations  are  solved  in 


V,  r(x,y) 


Schematic  structure  of  the  four-layer  and  two-layer  models 
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the  lowest  layer  of  the  four-layer  model,  values  of  the  variables  appropriate  for 
this  sigma  level  are  formed  by  /Y-weighted  vertical  averages  of  the  two-layer  model 
values.  These  averaged  values  are  the  ones  used  in  the  vertical  differencing  of 
the  four-layer  model  and  represent  the  means  of  feedback  from  the  two-layer  to  the 
four-layer  model. 

A  staggered  grid  is  used  in  both  the  vertical  and  horizontal  directions.  In 
the  vertical,  all  variables  are  layer  quantities  except  vertical  velocities,  which  are 
defined  at  interface  levels  (see  Fig.  1).  Horizontally,  velocities  are  defined  on 
staggered  points  (“x"  points  in  Fig.  21  which  surround  the  points  on  which  all  other 
variables  are  defined  ("o’’  points  in  Fig.  2).  In  order  to  increase  the  overall  model 
domain  size  and  move  the  lateral  boundaries  away  from  the  area  of  primary  interest, 
a  horizontal  nesting  of  the  model  is  employed  as  developed  by  Zhang  et  al.  (1986). 
The  horizontal  gridpoint  structure  of  the  nested  model  domains  is  shown  in  Fig.  2. 
A  fine  grid  mesh  (FGM)  with  20  km  resolution  is  nested  in  a  coarse  grid  mesh  tCGM) 
with  60  km  resolution.  A  3:1  ratio  of  FGM  points  to  CC-M  points  is  necessary  with 
a  staggered  grid  so  that  both  "x"  and  “o”  points  can  be  coincident  in  the  overlap 
region  'Zhang  et  al.  1986).  The  COM  domain  covers  1  320  km  x  1320  km  while  the 
FGM  domain  is  480  km  x  480  km  for  “o”  points  (all  displays  will  be  made  on  “o" 
point  arrays,  with  any  displayed  velocities  being  averaged  to  these  points).  The 
two-way  interactive  nesting  procedure  of  Zhang  et  al.  (1986)  is  used  with  a  few 
minor  modifications.  First,  in  the  calculation  of  tendencies  Tor  the  “o”  points  in  the 
FGM,  a  simple  linear  interpolation  is  used  between  CGM  points  nearest  the  boundary 
FGM  point  rather  than  the  “l.agrangian  interpolation"  used  by  Zhang  et  al.  (1986). 
Second,  no  additional  eddy  diffusion  is  used  in  the  region  of  the  boundary  at  this 
stage  of  model  development.  This  is  a  temporary  modification  so  that  any  noise 
resulting  from  improper  coding  or  inappropriate  nesting  criteria  will  be  more 
apparent.  Finally,  the  Newtonian  damping  scheme  applied  by  Zhang  et  al.  (1986)  near 
the  interface  is  not  yet  being  employed  for  the  same  reason.  We  fully  expect  that 
the  damping  scheme  and  perhaps  increased  diffusion  will  be  desirable  in  the  final 
model  to  help  reduce  noise  resulting  from  incompatabilities  between  the  CGM  and 
FGM  sol ut vins. 

The  model  with  nesting  as  described  above  is  referred  to  as  the  horizontally 
and  vertically  nested  model  (HVN).  Prior  to  completing  a  model  version  with  the 
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variables  are  defined  on  “o"  points. 
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two-way  interactive  nesting  of  Zhang  et  al.  (1986),  a  20  km  resolution  model  which 
employed  only  the  vertical  nesting  was  developed.  This  model,  which  had  a  domain 
size  of  520  km  x  520  km  (for  “o"  points)  is  referred  to  as  the  vertical  two-way 
interactive  nesting  model  (VTWIN).  Results  of  tests  with  both  models  will  be 
presented  in  section  3. 

It  is  possible  to  envision  another  degree  of  nesting  beyond  that  in  the  HVN 
model  which  results  in  a  considerable  reduction  of  grid  points,  and  hence  a 
significant  savings  in  computation  time.  Here,  we  retain  the  20  km  resolution  FGM 
shown  in  Fig.  2  only  for  the  two-layer  model.  The  four-layer  model  uses  a  60  km 
grid  spacing  over  the  entire  CGM  domain  including  the  region  "over"  the  FGM  two- 
layer  model.  Philosophically,  this  structure  represents  a  desire  to  achieve  high 
resolution  for  boundary  layer  processes  and  acknowledges  the  much  higher  density 
of  surface  observations  compared  to  upper  air  data.  The  vertical  two-way 
interactive  nesting  is  considerably  more  complicated  for  this  configuration  because 
it  requires  interpolation  of  coarse  four-layer  data  to  FGM  points  to  serve  as 
boundary  values  at  the  <JKt,m  level  and  subsequent  horizontal  (as  well  as  H-weighted 
vertical)  averaging  of  two-layer  variables  to  CGM  points  for  feedback  to  the  four- 
layer  model.  This  nesting  configuration  results  in  an  over-specification  of  the 
surface  pressure  variable,  r,  on  the  CGM  points  located  coincident  with  FGM 
points.  This  over-specification,  however,  is  no  different  from  that  which  is  present 
on  CGM  points  in  the  interface  region  when  the  horizontal  nesting  procedure  of 
Zhang  et  al.  (1986)  is  used  in  a  traditional  nested  model  (or  in  HVN).  Thus,  it  is 
hoped  that  with  careful  implementation  the  over-specification  will  be  slight  enough 
to  allow  successful  integration.  The  nesting  configuration  described  here,  and 
referred  to  as  horizontal  and  vertical  nesting  with  no  four-layer  FGM  (HYNNOFF), 
has  not  yet  been  tested.  Its  improvement  in  computation  speed  over  the  HVN  model 
would  be  the  result  of  a  substantial  decrease  in  gridpoints  and  also  because  CGM 
points  (and  hence  the  entire  four-layer  model)  are  only  integrated  once  for  every 
three  time  steps  of  the  FGM  points. 

r.  Other  numerical  details 

Time  integration  for  the  model  is  performed  using  the  leapfrog  scheme  with 
an  Asselin  filter.  The  time  step  for  FGM  points  is  20  s  and  for  CGM  points  is  60  s. 


!0 


The  flow  relaxation  condition  of  Davies  (1976)  is  used  on  lateral  boundaries 
following  the  work  of  Seitter  (1987)  who  found  that  this  condition  was  well-behaved, 
provided  a  simple  means  of  allowing  external  information  to  be  introduced  into  the 
model,  and  did  not  require  the  smoothing  operator  necessary  in  the  Perkey  and 
kreitzberg  (1976)  sponge.  The  flow  relaxation  condition  requires  a  5  gndpoint  wide 
region  near  the  boundary  for  application,  and  solutions  in  this  "relaxation  region” 
should  be  considered  modified.  For  the  non-nested  VTWIN  model,  this  means  that 
only  the  central  16  x  16  gridpoint  mesh  should  be  considered  as  having  a  true 
physical  solution.  For  the  HVN  model,  the  relaxation  region  is  contained  in  the 
CGM  collar  outside  the  CGM  points  providing  boundary  tendencies  to  the  FGM. 
Therefore,  the  entire  24  x  24  FGM  mesh  can  be  viewed  as  a  physical  solution. 

The  last  terms  on  the  RHS  of  (2.2),  (2.3)  and  (2.S)  include  a  "friction"  term, 
F,  which,  above  the  boundary  layer,  is  given  by  horizontal  eddy  diffusion.  This 

term  is  modeled  by  a  simple  Fickian  diffusion,  F'7‘0,  where  K  is  a  constant  eddy 

viscosity  and  0  is  the  variable  of  interest  (u,  v,  or  T).  In  most  simulations,  we  let 
K  —  2  x  10‘  m'  s~  for  the  momentum  components  and  K  =-  (1  for  temperature  (see 
section  3  for  a  more  complete  discussion  of  the  diffusion  of  temperature). 

Several  time-saving  approximations  are  being  introduced  into  the  coding  of 
the  model.  For  example,  energy  conservation  requires  that  the  values  of 

temperature  used  in  the  vertical  advection  term  of  the  thermodynamic  equation  be 
found  by  taking  the  mean  of  the  potential  temperature  (Anthes  and  Warner  1978). 
This  requires  the  calculation  of  potential  temperature  at  each  mid-layer  level, 

averaging  these  quantities  to  find  an  interface  mean  potential  temperature,  then 
using  the  Poisson  equation  to  calculate  the  resulting  interface  temperature. 
Normally,  this  requires  two  calculations  of  quantities  raised  to  non-integer  power 
for  each  model  layer,  for  each  grid  column,  per  time  step.  However,  a  small  bias  in 
the  pressure  level  determination  does  not  affect  the  potential  temperature 
calculation  substantially  and  that  same  bias  is  removed  when  conversion  back  to 
temperature  occurs.  Therefore,  it  is  possible  to  calculate  constant  factors  for  each 
level  at  model  initialization  based  on  average  surface  pressure  and  use  these  factors 
in  the  average  temperature  calculation  throughout  the  integration.  This  replaces 
the  non-integer  power  calculation  with  one  multiplication  while  producing  no 
stgnif icant  error. 
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The  VAX  VMS-Fortran  compiler  has  a  fairly  effective  optimization 
capability,  and  the  model  code  is  being  written  to  take  advantage  of  the  optimizing 
routine  as  much  as  possible.  This  includes  nesting  DO-loops  in  a  way  that  allows 
the  machine-language  code  to  employ  auto-increment  addressing  and  grouping 
variables  to  allow  some  precalculation  at  compilation  time  as  well  as  register  storage 
at  execution  time.  The  32-bit  word  length  of  the  VAX  has  shown  up  in  noticable 
roundoff  error  in  certain  calculations.  In  several  places  in  the  model,  special 
analytically  equivalent  forms  of  the  equations  have  been  used  to  minimize  the 
impact  of  the  roundoff  error  so  that  it  cannot  become  a  small  but  accumulating 
error  in  the  model.  One  place  where  this  showed  up  quite  noticeably  was  in  the 
application  of  the  flow  relaxation  condition  where  even  quiescent,  horizontally 
homogeneous  fields  were  being  modified  in  the  relaxation  region  unless  the  terms  in 
the  time-stepping  equations  were  written  to  reduce  their  sensitivity  to  roundoff 
error. 

3.  Tests  with  the  Model 

a.  Tests  using  VTWIN 

Extensive  tests  were  carried  out  with  the  non-nested  model  referred  to  as 
the  VTWIN  model.  For  all  the  tests  described  here,  the  model  had  no  radiation  or 
boundary  layer  parameterizations.  The  variable  boundary  layer  height,  a„,  was  set 
to  a  constant,  normally  halfway  between  the  surface  and  :t>„.  Inspection  of  the 
equations  given  in  section  2  reveals  that  the  specification  of  as  a  constant 
reduces  the  two-layer  model  to  a  regular  <7 -coordinate  model.  The  tests  described 
here,  therefore,  were  primarily  designed  to  serve  three  purposes:  1)  verify  the 
basic  coding  of  the  model;  2)  test  the  two-way  vertical  nesting  scheme  between  the 
two-layer  and  four-layer  models;  3)  assess  the  impact  of  the  coarse  vertical 
resolution  on  the  model’s  ability  to  simulate  flow  over  complex  terrain.  Some 
additional  tests,  not  described  below  in  detail,  were  made  with  o ^  set  so  that  one 
layer  of  the  two-layer  model  was  considerably  thicker  than  the  other.  These 
simulations  verified  the  consistency  of  the  H-weighted  vertical  averages  used  to 
produce  quantities  required  from  the  lowest  layer  of  the  four-layer  model  since 
changes  in  oh  did  not  adversely  affect  the  four-layer  model  solution. 
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After  a  series  of  geostrophic  adjustment  simulations  with  no  terrain  (similar 
to  those  of  Seitter  (1987)|,  a  mountain  ridge  was  added  and  tests  were  carried  out  on 
the  model's  ability  to  reproduce  mountain  lee-wavc  phenomena.  The  extremely 
coarse  vertical  resolution  in  the  upper  portion  of  the  model  was  expected  to 
severely  limit  its  ability  to  reproduce  mountain  lee-waves.  Interestingly,  many 
features  of  lee  waves  were  present  in  the  simulations  despite  the  inability  of  the 
model  to  properly  resolve  details  of  these  features.  Before  summarizing  the  results 
of  these  tests,  we  will  discuss  a  source  of  error  in  a-coordinate  models  which  does 
not  appear  to  have  an  acceptable  solution. 

During  the  lee-wave  testing,  it  became  clear  that  extreme  care  must  be 
exercised  when  applying  an  eddy  diffusion  to  variables  in  a  o-coordinate  model. 
Consider,  for  example,  the  eddy  diffusion  term  applied  to  temperature.  If  this  is 
handled  as  a  simple  Fickian  diffusion  in  the  form  KV‘T,  with  the  Laplacian 

evaluated  on  o-surfaces,  terrain  features  will  induce  an  accumulating  error  which 
may  eventually  destroy  the  simulation  through  dynamical  feedback.  To  see  how 
this  happens,  imagine  a  cr-surface  with  a  gridpoint  directly  above  a  peak.  All 

surrounding  gridpoints  on  the  same  model  surface  will  be  at  a  lower  physical  height 
because  the  cr-surface  will  deform  with  the  terrain  feature.  liven  if  the 
temperature  field  is  horizontally  homogeneous  (on  constant  height  surfaces),  a 
typical  lapse  rate  with  temperature  decreasing  with  height  will  produce  a  somewhat 
higher  value  of  T  on  the  gridpoints  surrounding  the  peak  compared  to  that  of  the 
peak  itself.  Application  of  the  Laplacian  will  produce  a  tendency  in  the 

thermodynamic  equation  which  will  raise  the  temperature  on  the  point  above  the 

peak  because  the  diffusion  term  will  tend  toward  no  gradients  on  sigma  surfaces. 
This  is,  of  course,  an  erroneous  tendency  which  is  an  artifact  of  the  simple 
treatment  of  the  diffusion  terms.  If  there  is  flow  over  the  peak,  which  will  tend  to 
cool  the  air  above  the  peak  adiabatically,  the  error  produced  by  the  diffusion  term 
will  accumulate  and  prevent  a  steady  state  from  being  reached.  This  problem  with 
the  eddy  diffusion  terms  will  be  most  noticeable  when  trying  to  simulate  nearly 
steady  flow  over  complex  terrain,  so  it  surfaced  early  in  test  simulations  of  lee- 
wave  phenomena. 

Ihe  best  way  to  reduce  errors  produced  by  the  diffusion  terms  is  not  clear. 
One  possibility  would  be  to  calculate  diffusion  on  a  constant  pressure  surface 
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through  the  gridpoint.  This,  however,  requires  an  interpolation  of  the  surrounding 
gridpoints'  data  from  their  c-surface  to  the  appropriate  pressure  surface  and  large 
truncation  errors  are  likely  to  result  —  especially  if  the  model  surfaces  are  not 
closely  spaced  in  the  vertical.  Mellor  and  Blumberg  (1985)  suggest  diffusion  on  cr 
surf  a  ces,  but  only  if  K  is  a  function  of  velocity  so  that  K  reduces  to  zero  for  zero 
velocity.  This  is  the  approach  used  in  some  other  mesoscale  models  (e.g.  Anthes 
and  Warner  1938).  ft  is  easy  to  see  that  this  does  not  eliminate  the  problem  of 
error  growth  when  the  flow  is  nearly  steady,  unless  the  atmospheric  lapse  rate  and 
the  slope  of  a  a  layer  combine  to  yield  a  constant  temperature  on  a  cr-surface.  The 
diffusion  term  could  be  reformulated  in  terms  of  potential  temperature.  This  would 
reduce  the  problem  in  some  types  of  flows  but  not  for  others,  and  would  require  a 
significant  increase  in  computation  for  these  terms.  We  anticipate  that  we  will 
continue  to  apply  a  simple  Fickian  diffusion  on  cr-surfaces,  but  that  we  will  set  the 
diffusion  coefficient,  K,  to  be  much  smaller  for  all  trace  variables  than  the  value 
used  for  momentum. 

Unlike  lee  wave  tests  performed  with  other  mesoscale  models  described  in  the 
literature  (Anthes  and  Warner  1978;  Nickerson  et  al.  1986),  the  tests  described  here 
were  performed  with  the  full  three-dimensional  model  rather  than  a  two-dimensional 
analog.  A  two-dimensional  analog  model  requires  a  different  gridpoint  staggering 
since  all  points  are  in  the  plane.  This  means  that  the  analog  model  is  an  inherently 
different  model  than  the  three-dimensional  one.  We  feel  it  is  better  to  test  the  full 
model  directly.  A  mountain  ridge  was  placed  in  the  terrain  field  running 
north-south  in  the  domain.  The  ridge  height  was  reduced  smoothly  to  zero  in  the 
vicinity  of  the  north  and  south  boundaries  to  reduce  inconsistencies  between  the 
steady  solution  over  the  ridge  and  the  boundary  values  which  are  held  constant  in 
the  flow  relaxation  condition.  This  resulted  in  a  ridge  that  was  independant  in  the 
y-direction  for  the  center  12  gridpoints  in  the  domain.  Solutions  were  displayed  on 
an  east  west  cross  section  through  the  center  of  the  domain.  The  resulting  flow 
should  be  similar  to  two-dimensional  solutions  in  the  literature  (though  not 
suffering  from  unrealistic  two-dimensional  constraints).  These  simulations  provided 
a  stringent  test  for  the  model. 

Figure  3a  shows  the  potential  temperature  field  after  1 2  h  of  simulation  for  a 
1  km  Gaussian  profile  ridge.  The  u-component  velocity  field  for  this  simulation  is 
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FIG.  3.  Fast-west  vertical  cross-section  through  the  center  of  the  VTWIN  model  at 
13  h  simulated  time  showing  (a)  potential  temperature  and  (b)  u-component 
velocity.  The  scale  along  the  left  shows  height  in  km,  lsentropes  in  panel  (a) 
are  labeled  in  K  along  the  right  boundary. 
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shown  in  Fig.  3b.  1  he  l-.S.  Standard  Atmosphere  temperature  profile  and  20  m  s_1 

/^-velocity  at  all  levels  served  as  the  initial  conditions  for  these  simulations  and 
these  values  are  held  constant  on  all  boundaries  by  the  flow  relaxation  scheme. 
The  Coriolis  parameter  was  set  to  zero  for  these  simulations  in  order  to  allow 
comparison  between  these  tests  and  other  published  results.  In  Fig.  3a,  the  points 
show  the  positions  of  model  gridpoints,  with  the  plotted  lsentropes  resulting  from 
linear  interpolation.  In  Fig.  3b,  the  velocity  values  are  plotted  at  each  gridpoint 
location  and  ttie  coutours  have  been  subjectively  analyzed.  Clearly  the  coarse 
vertical  resolution  leads  to  a  smoothing  of  the  lee-wave  structure,  but  comparison 
with  a  similar  simulation  by  Nickerson  et  ai.  (198b,  his  Fig.  4a)  shows  that  the  major 
features  are  present  especially  in  the  lower  6  km  where  the  model  resolution  is 
better.  A  velocity  maximum  of  over  30  ms  is  located  near  the  surface  just 
downwind  of  the  peak  and  a  velocity  minimum  of  less  than  8  m  s  ’  is  located  at 
midlevels  directly  above  the  low  level  maximum.  The  velocity  structure  and 
potential  temperature  fields  both  show  a  distinct  upwind  tilt  for  the  waxes.  There 
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FIG.  4.  Average  rate  of  change  of  surface  pressure  for  12  h  simulation  shown  in 
Fig.  3. 
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is  no  damping  layer  in  the  upper  portion  of  the  model,  so  the  amplitude  of  the  wave 
is  larger  there  than  it  should  be.  Unlike  the  Nickerson  et  al.  (1986)  solution,  the 
model  has  no  frictional  boundary  layer  for  these  tests,  so  the  velocities  are  carried 
down  to  the  surface. 

Figure  4  shows  a  plot  of  the  average  rate  of  change  of  the  surface  pressure 
for  the  1_'  h  of  the  simulation  shown  in  Fig.  3.  The  dashed  curve  is  the  average  of 
■  cfT  dt  over  the  domain  before  the  application  of  the  flow  relaxation  condition.  A 
zero  value  for  this  quantity  indicates  perfect  dynamic  balance  within  the  model 
domain.  It  does  not  reach  perfect  balance  because  the  boundary  values  are  not 
exactly  consistent  with  the  lee-wave  solution.  The  solid  curve  represents  the 
average  over  the  domain  of  . d -r  dt 1  after  the  application  of  the  boundary  condition 
(see  Seitter  (1987)  for  a  more  complete  discussion  of  the  difference  between  these 
two  quantities).  A  value  of  zero  for  this  quantity  indicates  a  perfectly  steady 
solution.  This  figure  indicates  a  period  of  adjustment  within  the  model  for  the 
first  hour  during  which  the  flow  becomes  nearly  balanced  (except  for  the  gridpoints 
in  the  flow  relaxation  region  near  the  boundaries).  After  the  first  hour  the  model 
solution  is  very  steady. 

Other  test  simulations  with  lee-wave  phenomena  included  varying  the  width 
and  sharpness  of  the  mountain  ridge,  varying  the  stability,  and  varying  the  eddy 
coefficient.  The  solution  changed  in  only  minor  ways  with  changes  in  the  mountain 
shape  or  atmospheric  stability  and  always  in  ways  consistent  with  lee-wave  theory. 
Reduction  of  the  eddy  diffusion  coefficient  for  momentum  by  a  factor  of  10  or 
more  resulted  in  short  wavelength  noise  being  generated,  especially  upwind  of  the 
mountain.  This  is  thought  to  be  primarily  an  artifact  of  the  closeness  of  the 
upwind  boundary  and  interactions  between  the  boundary  condition  and  the  flow. 
However,  Anthes  and  Warner  (1978)  indicate  that  low  resolution  in  the  vertical 
requires  a  nondimensional  eddy  coefficient  of  about  0.01  in  order  for  a  reasonable 
solution  to  be  produced  when  no  upper  damping  iayer  is  present  in  the  model.  The 
solution  shown  in  Fig.  3  was  produced  with  a  nondimensional  eddy  coefficient  of 
0.005. 

The  solutions  from  these  tests  are  quite  good  in  the  lower  troposphere  where 
the  emphasis  of  the  model  is  being  placed,  nearly,  the  model  is  incapable  of  being 
used  to  investigate  the  details  <f  lee-wavc  phenomena,  but  we  are  very  encouraged 


MCI.  5.  last-west  vertical  cross-section  through  the  center  of  the  HYN  model  fine 
grid  mesh  at  12  h  simulated  time  showing  (a)  potential  temperature  and  (b)  u- 
component  velocity. 
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by  its  ability  to  capture  those  features  which  will  impact  the  forecast  at  lower 
levels  during  these  events. 

b.  Test s  with  the  H t-'.V  model 

Mountain  lee-wave  simulations  similar  to  the  ones  described  above  (as  well  as 
geostrophic  adjustment  and  other  tests)  have  also  been  performed  with  the 
horizontally  nested  HYN  model.  As  in  the  non-nested  simulations,  a  north  south 
mountain  ridge  of  1  km  height  was  placed  in  the  domain.  The  northern  and 
southern  ends  of  the  ridge  sloped  smoothly  to  zero  height  and  the  ridge  was 
designed  to  fit  in  only  the  PGM  domain  (though  its  northern  and  southern  ends 
extended  somewhat  into  the  overlapping  feedback  region  of  the  CGM).  The 
resulting  ridge  had  a  constant  height  over  10  gndpoints  in  the  y-direction. 
Simulations  out  to  1 2  h  were  carried  out  for  several  ridge  profiles. 

Figure  5a  shows  an  east  west  cross-section  of  potential  temperature  at  12  h 
in  the  FGM  domain  for  the  same  Gaussian  profile  ridge  as  shown  in  Fig.  3a.  These 
results  appear  very  similar  to  those  of  the  non-nested  model.  Despite  the  lack  of 
diffusion  on  temperature  and  no  Newtonian  damping  term  near  the  boundary,  the 
solution  is  quite  well-behaved.  The  cross-section  of  u-velocity  for  this  simulation, 
shown  in  Fig.  5b,  shows  that  the  lee-wave  has  a  smaller  amplitude  than  that  of  the 
the  non-nested  simulation.  This  appears  to  be  related  to  the  close  lateral 
boundaries  in  the  non-nested  model  which  held  the  velocity  fixed  at  JO  m  s  ‘.  The 
nested  model  allows  the  velocity  to  adjust  at  the  FGM  boundary  and  it  appears  that 
the  wave  took  on  a  broader  structure  of  lower  amplitude  when  not  constrained  by 
the  flow  relaxation  condition  on  FGM  boundaries. 

Fo  show  the  smoothness  of  the  solution  at  the  interface  between  the  CGM 
and  FGM,  Fig.  6  shows  a  vertical  cross-section  of  potential  temperature  through  the 
CGM  domain.  The  extent  of  the  FGM  domain  is  shown  by  the  vertical  dashed  lines, 
and  the  values  plotted  within  this  region  were  obtained  by  applying  the  same  9-point 
operator  that  is  used  to  provide  feedback  from  the  FGM  to  the  CGM  in  the  two- 
way  interactive  nested  procedure  (Zhang  et  al.  1986).  Comparison  of  Figs.  5a  and  6 
show  that  the  lee-wave  structure  would  not  be  resolved  very  well  by  a  60  km 
resolution  model.  This  did  not,  however,  lead  to  an  incompatibility  between  the 
FGM  and  CGM  solutions  at  the  interface,  so  no  discontinuities  are  apparent. 
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FIG.  6.  East-west  vertical  cross-section  though  the  center  of  the  HVN  model  coarse 
grid  mesh  at  12  h  simulated  time  showing  potential  temperature. 

Figure  7  shows  the  variation  with  time  of  the  FGM  domain  average  rate  of 
change  of  surface  pressure,  idir  dl  ,  for  the  12  h  simulation  shown  in  Fig.  5. 
Comparison  with  Fig.  4  shows  that  the  nested  model  requires  a  somewhat  longer  time 
to  reach  a  nearly  steady  state  (about  2  h)  and  undergoes  larger  oscillations  while 
adjusting.  This  is  clearly  related  to  the  lack  of  damping  in  the  FGM  of  the  nested 
simulation  compared  to  that  provided  by  the  flow  relaxation  condition  in  the  non¬ 
nested  model.  It  would  also  appear  that  Lamb  waves  produced  during  the  adjustment 
period  suffer  partial  reflection  at  the  FGM-CGM  interface  due  either  to  aliasing  of 
the  waves  or  as  a  result  of  the  over-specification  of  the  pressure  in  the 
FGM-CGM  overlap  region.  Still,  the  model  appears  quite  stable  —  even  this 
version  which  has  little  means  of  removing  noise  in  the  FGM.  This  is  very 
encouraging  and  suggests  that  even  though  the  coarse  vertical  resolution  of  the 
model  fails  to  adequately  resolve  the  details  of  the  lee-wave  phenomena,  there  is  no 
reason  to  fear  that  topographically  induced  waves  will  lead  to  disturbances  which 
will  destroy  other  aspects  of  the  simulation. 
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TIME  (hr) 

FIC.  7.  Average  rate  of  change  of  surface  pressure  in  PGM  for  12  h  simulation 
shown  in  Fig.  5. 

Despite  the  considerable  increase  in  the  size  of  the  physical  domain  for  the 
nested  model  and  the  added  complexity  of  the  horizontal  nesting  procedure,  the  HVN 
model  runs  in  nearly  the  same  amount  of  time  as  the  non-nested  VTWIN  model  (just 
under  1  h  CPU  time  for  a  12  h  simulation  on  a  MicroVAX  3).  This  is  partly  a 
result  of  coding  improvements  made  to  the  HVN  model  which  were  not  utilized  in 
VTWIN.  Mostly,  however,  the  comparable  speed  is  a  result  of  two  apsects  of  the 
nesting  itself.  First,  the  FGM  of  the  HVN  has  only  a  24  x  24  gridpoint  domain  (for 
"o"  points)  compared  to  the  26  x  26  domain  of  the  VTWIN'  model,  so  HVN  has  about 
15%  fewer  gridpoints  at  20  km  resolution.  Second,  the  60  km  resolution  CGM  points 
surrounding  the  FGM  in  HVN  are  only  integrated  once  for  every  3  timesteps  of  the 
FGM  (and  the  CGM  collar  only  requires  integration  on  about  2/3  as  many  points  as 
the  FGM).  It  seems  clear  that  the  nested  HVN  model  is  preferable  to  the  non¬ 
nested  VTWIN,  and  does  not  increase  computation  time  enough  to  threaten  the 
operational  goal  of  this  project. 
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4.  Boundary  layer  parameterization 

All  of  the  test  simulations  described  in  the  previous  section  were  carried  out 
in  versions  of  the  model  which  had  a  constant  boundary  layer  height,  cr.,.  As  stated 
in  section  2,  however,  a  major  component  of  this  work  is  the  development  of  a  model 
with  variable-height  boundary  layer.  The  formulation  of  the  boundary  layer 
parameterizations  which  will  be  incorporated  into  the  model  described  above  is 
presented  in  this  section.  The  boundary  layer  parameterization  is  composed  of 
several  parts:  radiation;  surface  energy  and  moisture  balance;  boundary  layer 

height;  and  boundary  layer  air  temperature  and  specific  humidity.  Each  of  the 
formulations,  described  in  the  following  subsections,  is  designed  to  balance 
computational  speed  with  physical  accuracy.  In  most  cases,  the  physical  accuracy 
is  restricted  by  the  limited  vertical  resolution  of  the  model  rather  than  by  the 
mathematical  formulation. 

u.  rcidi'ition 

The  radiation  parameterization  is  taken  from  katayama  (1972)  and  is  a  routine 
originally  designed  for  use  in  the  UCLA  GCM.  The  incident  radiation  and  infra-red 
(IR)  emission  are  calculated  separately.  The  model  incorporates  an  exponential  fit 
to  the  data  for  specific  humidity  to  allow  simple  integration  of  water  content.  CO: 
is  included  in  a  fixed  form  based  on  experimental  data,  and  its  contribution  is  then  a 
constant. 

The  influx  of  radiation  is  computed  by  starting  with  the  solar  constant  and 
modifying  it  for  albedo  at  the  top  of  the  atmosphere.  Scattered  and  absorbable 
radiation  are  computed  separately,  the  fraction  being  assumed  constant  (35% 
available  for  absorption,  65%  scattered  to  the  ground).  T  he  scattered  part  of  the 
incident  radiation  is  corrected  for  multiple  reflection  between  the  atmosphere  and 
the  ground  and  given  by 

CLW s  -  ( .651  ).S'(jcos  Z'/' i  1  I  (4.1) 

m 

where 

S0  =  solar  constant  as  a  function  of  day  of  year 


ZT  —  zenith  angle  for  time  of  day  and  location 

as  =  scattering  albedo  for  the  atmosphere  (if  clouds  are  present  they 
determine  the  scattering  albedoi 

a,.  —  albedo  of  ground  surface  as  a  function  of  hour  angle  and  surface 
characteristics  (Wetzel  l^S) 

If  a  cloud  layer  is  present,  its  presence  is  felt  by  both  scattered  and  absorbable 
components.  If  the  cloud  is  thick  enough,  and  covers  enough  sky,  incident  radiation 
can  be  shut  off.  The  model  allows  for  variable  amounts  of  cloud  in  each 
atmospheric  layer  expressed  as  a  percentage.  Only  one  layer  of  cloud  is  allowed, 
but  it  may  be  composed  of  one  or  several  atmospheric  layers  of  various  percentage 
coverage. 

Fractional  absorption  by  water  vapor  is  calculated  by 

E  ,U.303 

ABSi  (t).rn  h-%L- 
cosiZT 


(4.2) 


where  EH  ,  is  the  effective  amount  of  water  vapor  in  layer  i.  The  radiation  which 
finally  is  absorbed  in  the  soil  becomes  one  component  of  the  surface  energy  balance. 
The  absorbed  part  of  the  incident  radiation  at  the  ground  is 


GLW  i  =-=-  lO.349  Sccos.ZT'  -  ^ABSti. 
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(4.3) 


where  the  sum  is  taken  over  all  layers.  The  total  absorption  at  the  ground  is  then 


GAB  —  il  -  cijnCl.W * 


GLW 


(4.4) 


To  find  the  IR  flux,  the  equation  of  radiative  transfer  is  solved  subject  to 
the  boundary  conditions  that  downward  IR  flux  at  the  top  of  the  atmosphere  is 
zero,  and  the  upward  IR  flux  at  the  earth’s  surface  is  the  black-body  radiation  at 
the  surface  temperature.  Weighted  transmission  functions  are  used,  corrected  for 
the  pressure  dependence  of  absorption  by  defining  an  effective  amount  of  an 
absorber.  The  total  transmission  function  is  assumed  to  be  the  product  of  the 
individual  ones  for  CO  and  H  0.  Downward  flux  is 
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where 

-rB,  —  <T7'f 

cr  =  Stephan-Boltzmann  constant 

r,  t  =r  mean  total  trnasmission  functions  for  effective  absorber  u  at 
temperature  T 

TV  —  critical  temperature  which  divides  the  region  of  weak  temperature 
dependence  of  r  to  that  of  strong  dependence  of  t. 


The  weak  region  is  210-320  K  for  water  vapor,  so  letting  T-  —  220  K,  the  weak 
dependence  region  need  only  have  a  mean  temperature  specificied  (T).  Similarly,  the 
upward  flux  is 


I  Ri. 
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(4.6) 


and  the  net  upward  flux  is 


IR i  —  IRu  --  IRj 


(4.7) 


The  only  difficulty  is  determining  the  proper  transmission  function  near  the 
particular  level  where  r  varies  exponentially.  The  model  uses  an  interpolation 
factor  which  is  an  empirical  function  of  pressure,  mixing  ratio  and  layer  thickness. 
This  allows  proper  calculation  of  T  without  a  fine  vertical  mesh.  The  mean 
transmission  functions  are  defined  by  empirical  formulae  at  T-  —  220  K  and  T  — 
260  K.  Temperature  dependence  of  'for  CO  is  neglected,  so  a  mean  T  for  CO:  is 
used  based  on  pressure  and  amount  of  CO,.  The  distribution  of  CO>  at  each  level 
is  a  constant. 

The  IR  flux  is  computed  only  for  the  surface,  since  the  IR  cooling  rates  in 
the  free  atmosphere  are  insignificant  on  diurnal  time  scales.  This  saves 
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considerable  computation  time.  The  net  radiation  is  then 

SR  =  GAB  IR j  .  (4.8) 

Radiation  calculations  were  made  with  the  above  scheme  for  a  varying 
number  of  atmospheric  layers.  These  calculations  showed  that  when  layers  were 
thicker  than  100  mb,  errors  occurred  in  the  net  radiation  values.  Comparisons  were 
made,  for  example,  using  a  sounding  that  originally  had  19  unevenly  spaced  levels. 
Reducing  this  to  only  7  levels  produced  NR  values  in  error  by  nearly  15%.  If, 
however,  the  routine  started  with  the  same  7  levels  but  was  allowed  to  linearly 
interpolate  a  new  level  in  the  middle  of  any  layer  thicker  than  100  mb,  the  error  in 
SR  was  reduced  to  less  than  10%  —  despite  the  fact  that  no  additional  vertical 
resolution  in  temperature  was  available. 

b.  surface  energy  balance 

The  surface  energy  balance  has  the  form 

NR  =  SH  +  LH  +  GS  (4.9) 

where  NR  is  the  net  radiation  incident  on  the  surface  (as  defined  in  the  previous 
subsection),  SH  is  the  sensible  heat  flux  upward  from  the  surface,  I.H  is  the  latent 
heat  flux  upward  from  the  surface,  and  GS  is  the  soil  heat  flux  downward  into  the 
ground  (heating  the  soil).  The  NR  value  is  determined  from  the  radiation 
parameterization  (see  eq.  (4.8)1.  The  other  terms  are  parameterized  as  follows. 

The  sensible  heat  flux  and  latent  heat  flux  (SH,  LH)  are  parameterized  using 
Monin-Obukhov  similarity  theory  for  the  planetary  boundary  layer  (PBL).  The 
fluxes  depend  on  the  vertical  gradients  of  temperature  and  specific  humidity  in  the 
surface  layer,  the  surface  geostrophic  wind,  and  the  stability  of  the  boundary 
layer.  The  theory  assumes  that  the  structure  of  temperature  and  moisture  in  the 
PB1  have  forms  which  can  be  described  by  universal  structure  functions  when 
scaled  equations  are  used.  (here  are  actually  two  structures  involved,  since  the 
PBL  contains  at  least  two  distinct  layers  the  surface  layer  and  the  boundary 
layer.  Although  the  two-layer  model  structure  does  not  include  an  explicit  surface 


layer,  one  is  assumed  to  be  present  by  the  parameterization.  This  “surface  layer" 
is  assumed  to  be  a  constant  5  mb  thick.  If  the  functions  are  required  to  be  matched 
at  their  common  boundary,  the  following  results  are  obtained 
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where 

Ugr,  —  components  of  the  surface  geostrophic  wind 
-  potential  temperature  at  the  top  of  the  PBL 
0  —  potential  temperature  at  the  ground 

0*  -  •  SH  u* 

um  •---  friction  velocity 
k  —  von  Karman’s  constant 
C  —  dimensionless  constant 
U  •=  kuM/fl. 
v,  -  h  /. 

/  —  Coriolis  parameter 
/.  -  Obukhov's  length 

z.  roughness  length  (function  of  location) 
h  -=  height  of  the  boundary  layer 


(4.14) 

(4.15) 
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A  Ci  ,u.  i  -  universal  functions  for  a  stable  boundary  layer 

—  universal  functions  for  an  unstable  boundary  layer 


where  the  universal  functions  are  those  given  by  Arya  (1975). 

To  determine  the  fluxes,  we  use  inverted  forms  of  these  equations  which  are 
based  on  two  parameters: 

stable  case  —  S  —  (4. IS) 
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unstable  case  — 
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where 

Go  —  surface  geostrophic  windspeed 
y  =  acceleration  of  gravity 
6  —  mean  potential  temperature. 

The  inverted  equations  have  the  form 
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(4.16) 
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where  r0  is  the  surface  stress  and  o  is  the  density  of  air.  The  latent  heat  flux 
( LH )  is  computed  by  assuming  that  it  obeys  the  universal  function  for  SH,  that  is. 


0-  -  0. , 

<7*  '  0* 

where  q  is  the  specific  humidity  and 


(4.21) 
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c.  ground  variables 

Ground  temperature  ( TG )  and  ground  wetness  (G’lV)  are  parameterized  by 
“force  restore"  methods  from  Bhumralker  (1975)  and  Deardorff  (1977).  Following 
Bhumralker  (1975),  heat  conduction  in  the  soil  is  described  by 

oT ,/  z,t  k  3Tf(z,lj  (4  ^3) 

r  iz: 


where 

Tg<z,ti  =  soil  temperature  at  depth  z,  time  t 
K  =  thermal  conductivity  of  soil 
c  =  volumetric  heat  capacity. 

We  assume  that  TG  is  described  by 

TG  —  T  *-  AT(,sin(wf '  (4.24) 

where 

T  =  average  temperature  of  the  soil,  assumed  to  be  invariant  with 
depth  (in  practice,  the  50  cm  soil  temperature  is  used) 
i-.T ,  =  amplitude  of  the  variance 
jj  —  frequency  of  the  variance. 
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Then  the  solution  of  (4.21)  is 

Tj.z,/1  =  T  -  \T  expi  -- z /d)[sin'u)/  z /d'1 

where  d  —  ‘2k  cu)'‘  *  is  the  depth  at  which  the  amplitude  of  2.  7' 
an  infinitely  thin  soil  layer,  the  heat  flux  into  the  soil  at  depth  3 
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.37Vz.li 

?z 


Combining  (4.25)  and  (4.2b)  gives 


G z,f  =  J^y^j  .' To  r  z  ^(sinUvl  z  d;  ♦  cosia >t-Z't 


Eliminating  2iT,,  one  obtains 
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Consider  a  layer  of  soil  from  the  surface  (z  --  0)  to  some 
rate  of  temperature  change  for  this  layer  is  given  by 
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If  the  approximation  is  made  that 


Ty-'Z=1  cm,  ti  i  TC 


then  (4.29)  becomes  (with  the  use  of  (4.9)1, 
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is  negligible.  For 
:  is 

(4.26) 

f  ti  .  (4.27) 

(4.28) 

depth  z.  The  time 

(4.29) 

(4.30) 

(4.31) 
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The  surface  soil  moisture  (GW)  is  found  by  assuming  that  it  changes  due  to 
three  main  processes  —  precipitation,  evaporation,  and  flux  from  below.  The  bulk 
soil  moisture  (CUB)  is  assumed  to  be  constant  over  the  period.  According  to 
Deardorff  (1977)  the  bulk  soil  moisture  changes  over  a  time  scale  of  a  few  weeks,  so 
GW  B  can  certainly  be  assumed  constant  for  a  24  h  period  with  little  loss  of 
accuracy.  The  surface  soil  moisture  is  changed  according  to 

3(;U  _  coLH  \  -  Pr )  c.tGWj—Clv  B)  (4  33) 
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where 

GW B  =-  percent  bulk  soil  saturation  (top  50  cm) 

GW  —  percent  surface  soil  saturation 
d,  -  depth  of  diurnal  cycle  (=  10  cm) 
k  =  latent  heat  of  evaporation 
p =  density  of  water  (1  gm  cm  ’) 

U'iwx  =  field  capacity  soil  moisture 
T  =  period  of  cycle 
c;,  c.  '■  nondimensional  constants 
P  —  precipitation  rate. 

Deardorffs  values  for  c,  and  c.  were  computed  from  data  of  Jackson  (1973), 
measurements  taken  over  bare  soil  near  Pheonix,  Arizona  in  March.  This  gives 
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(4.34) 


Notice  that  the  middle  value  of  c  is  a  linear  interpolation  between  the  two 


extremes. 
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d.  boundary  layer  height 


1)  unstable  boundary  layer 

The  unstable  PBL  is  assumed  to  be  well  mixed  below  an  inversion 
characterized  by  a  jump  discontinuity  in  potential  temperature  (TO).  The  depth  of 
the  unstable  boundary  layer,  h,  and  the  strength  of  the  inversion,  A9,  are  predicted 
according  to  Zeman  and  Tennekes  (1977).  Their  method  assumes  that  the  PBL  depth 
changes  due  to  turbulent  entrainment  of  air  above  the  inversion  into  the  PBL.  I  he 
energy  comes  from  the  virtual  SH  flux  at  the  surface,  and  the  change  of  depth 
with  time  depends  on  the  strength  of  the  inversion.  They  use  the  turbulent  kinetic 
energy  (TKE)  budget  to  develop  a  simple  set  of  equations  to  describe  this  process. 

The  sensible  heat  flux  at  the  inversion  is  equal  to  the  temperature  jump,  Lb, 
times  the  rate  of  rise  of  the  inversion 


l 'SH,  -  (4.35) 

where  V SH,,  is  the  virtual  sensible  heat  flux  at  the  inversion.  The  inversion 
strength  changes  as  a  function  of  entrainment  of  stable  air  from  above,  and  net 
sensible  heat  transfer  inside  the  boundary  layer.  It  is  give^by 


3J.0  y  dh 

a  t  1  a  t 


\  SH  LSH,  ‘h 


(4.36) 


where  L SH  is  the  virtual  sensible  heat  flux  at  the  surface  and  is  the  potential 
temperature  lapse  rate  above  the  PBL.  The  TKE  budget  can  be  written  as 
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which  can  be  expanded  into 
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where 
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wm  =  convective  velocity  scale  =  1  SH  hj 

cvbv  —  Brunt-Vaisala  frequency  =  J^j* 

TS  —  temperature  at  the  top  of  the  surface  layer 

and  Cf,  C f,  and  Cd  are  dimensionless  constants.  Substituting  for  9h/3t  from  (4.35) 
yields 

(4.39) 

Substituting  for  u>*  gives 
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(4.40) 


The  values  of  the  dimensionless  coefficients  CJt  Ct,  and  Ct  are  taken  to  be  (Zeman 
1975) 


Cd  =  0.024 

Cf  =  0.50  (4.41) 

C,  -  3.55 


We  can  write  the  rate  of  change  of  boundary  layer  height  as 
dh  I'SH* 

a?  ~  ‘  7J  •  (4-4^ 

In  the  case  where  A 9  ™  0,  no  inversion  exists  and  the  atmosphere  presents  no 
barrier  to  inversion  rise.  In  this  case,  the  model  assumes  a  very  small  value  for 
Afl,  since  the  inversion  must  rise  at  a  rapid  but  finite  rate  due  to  turbulent 
entrainment. 

2)  stable  boundary  layer 


The  depth  of  the  stable  boundary  layer  is  calculated  using  a 
parameterization  from  Yamada  (1979)  that  starts  with  the  thermal  energy  equation 
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for  flat  terrain, 
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(4.43) 


where  SHZ  is  the  sensible  heat  flux  at  level  z  and  (30  '3 Or  is  the  rate  of  change  of 
temperature  due  to  longwave  radiation.  Integrating  this  between  the  surface  and 
the  top  of  the  boundary  layer,  h,  yields 

f‘  f,  ^  -  SH,  -  SH.  +  fh  (|f)r  dz  (4.44) 

where  the  subscript  s  refers  to  the  surface  (actually,  the  top  of  the  surface  layer) 
and  subscript  h  refers  to  the  top  of  the  boundary  layer.  If  we  assume  that  the 
vertical  profile  of  potential  temperature  has  a  simple  cubic  form,  increasing  upwards 
(supported  by  observations),  the  equation  can  be  integrated  to  yield 


1  a  a  .dh 

4  05  cU 


f 


hi  :>3<h? 

4*  Of 


-I 


30,) 
dt  > 


(4.45) 


Since  turbulence  in  a  stable  PBL  decreases  rapidly  with  height,  at  the  top  of 
the  PBL  SH  —0,  so  (4.43)  becomes 


(30)  =  B0*. 

^3Mr  dt 


(4.46) 


At  the  surface,  the  cooling  of  the  air  is  taking  place  primarily  by  radiation,  so 
(4.43)  becomes 
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Linearly  interpolating  between  these  two  expressions  gives 
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Equation  (4.44)  then  becomes 


dz  =  SHs 
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(4.49) 


The  change  of  potential  temperature  above  the  PBL  is  much  smaller  than  the  change 
at  the  surface,  so  t90*/3f,  0.  Further,  SH*  is  typically  much  smaller  then  SH, 
it,  too,  can  be  set  to  zero  as  a  good  approximation.  Equating  (4.45)  and  (4.49),  we 
obtain  the  rate  equation  for  the  stable  PBL  depth 


0„  -0S  9f 
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e.  boundary  layer  variables 

The  surface  temperature  is  found  as  a  by-product  of  the  PBL  depth 
calculations.  For  the  unstable  PBL,  we  proceed  as  follows: 


L9p  =0,-4  {Qr  - 0 ,)  -  0r 


(4.51) 


where 

6P  ■=  potential  temperature  above  jump  discontinuity 
0,  —  potential  temperature  of  the  free  atmosphere  above  the  PBL 
zt  =  thickness  of  the  next  model  layer  above  the  PBL 
Lb  =  change  in  PBL  depth  over  one  time  step. 

The  potential  temperature  at  the  top  of  the  5  mb  surface  layer  is  given  by 


0s  Op  —  L9 


(4.52) 


which  allows  the  surface  temperature  to  be  found  by  Poisson’s  equation  as 


TS  -  0 
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(4.53) 
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For  the  stable  PBL,  we  start  in  a  similar  manner,  except  that  there  is  no 
inversion  jump  discontinuity.  We  simply  interpolate  in  potential  temperature 
between  the  PBL  top  and  the  ground  to  find  0S  then  compute  TS  with  (4.53). 

For  humidity,  a  moisture  budget  is  used  for  both  stable  and  unstable  PBLs. 
This  budget  can  be  written 

[flux  of  moisture!  ,  [entrainment  of  moisture!  U  ^41 

[  into  PBL  j  j  from  above  PBL  J  ‘ 

The  flux  of  moisture  is  given  by  LH,  and  the  entrainment  of  moisture  is  a  result  of 
the  growth  of  the  PBL  and  thus  a  function  of  Ah.  The  rate  of  change  of  specific 
humidity  given  by  (4.54)  is  then  used  in  the  model  equation  governing  the  change  of 
specific  humidity. 

5.  Conclusions  and  Future  Work 

The  problem  of  developing  a  rneso-#  numerical  model  capable  of  providing 
operational  forecast  guidance  while  running  on  a  super-micro  class  computer  is  not 
an  easy  one.  The  model  described  in  this  report  appears  to  be  one  possible  solution 
for  this  problem.  While  much  simpler  than  most  mesoscale  models,  it  still  attempts 
to  treat  important  physical  processes  in  a  realistic  way  —  though  much  differently 
from  other  mesoscale  models.  The  coarse  vertical  resolution  is  clearly  a  sacrifice, 
but  tests  indicate  that  the  model  is  stable  and  well-behaved  even  while  simulating 
flow  over  complex  terrain.  It  is  noteworthy,  in  regard  to  this  point,  that  early 
versions  of  the  PSU/NCAR  mesoscale  model  had  nearly  the  same  vertical  resolution 
as  the  model  developed  here. 

A  major  aspect  of  the  model  development,  the  variable  depth  boundary  layer 
formulation  of  the  two-layer  model,  has  not  been  tested  within  the  three-dimensional 
model.  One-dimensional  testing  of  the  boundary  layer  parameterizations  has 
produced  realistic  results,  but  examination  of  the  usefulness  of  the  ^-coordinate  as 
a  means  of  representing  a  variable-depth  model  layer  within  a  three-dimensional 
model  awaits  the  completion  of  the  integration  of  the  parameterization  into  the 
complete  model.  This  is  being  worked  on  at  the  writing  of  this  report,  and 
represents  the  first  stage  of  future  work  on  the  model. 
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Other  future  work  includes  the  addition  of  moisture  to  the  model.  The 
*  equations  for  humidity  are  being  added  at  the  same  time  as  the  boundary  layer  and 

radiation  parameterizations  so  that  the  humidity  variables  required  for  these 
a  packages  are  present.  Prognostic  equations  for  clouds  and  precipitation  will  not  be 

included  until  extensive  testing  for  clear-sky  conditions  has  been  completed.  These 
tests  will  include  such  phenomena  as  the  development  of  the  sea-breeze,  diurnal 
variation  of  the  boundary  layer,  and  diurnal  mountain-valley  circulations. 
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